Aplicación del análisis de componentes principales a las imágenes
Juan David Ospina Arango
Analítica Predictiva
Universidad Nacional de Colombia
2 de julio de 2020
Componentes principales y análisis de imágenes
En este documento se explora la aplicación del análisis de componentes principales sobre imágenes.
Las imágenes han sido obtenidas de A database of facial expressions in younger, middle-aged, and older women and men y solo se han descargado las 72 imágenes que se pueden descargar sin haberse registrado.
Para citar el uso de las imágenes se usa la siguiente referencia
Holland, C. A. C., Ebner, N. C., Lin, T., & Samanez-Larkin, G. R. (2019). Emotion identification across adulthood using the Dynamic FACES database of emotional expressions in younger, middle aged, and older adults. Cognition and Emotion, 33, 245-257. doi:10.1080/02699931.2018.1445981.
Carga de las imágenes
Se utilizarán algunas funciones de la librería imager para la carga y la transformación de las imágenes.
A continuación se carga la librería imager y se indica el directorio donde se encuentran las imágenes:
library(imager)
filenames <- list.files(path = "./img", pattern="*.jpg")El objeto filenames es una lista que contiene las imágenes. Para cargar una imagen con la funcón load.image() se debe indicar la ruta completa. Para completar el nombre de la ruta se usa la función paste0():
im <- load.image(paste0("./img/",filenames[1]))
plot(im)
Ahora se muestra la imagen en escala de grises con la función grayscale:
plot(grayscale(im))
La dimensión de la imagen es:
dim(im)## [1] 2835 3543 1 3Esto quiere decir que el objeto im es una sola imagen de \(2835 \times 3543\) pixeles, con tres canales de color (RGB).
Ahora se utiliza la función resize() para reducir la imagen a diferentes tamaños que corresponden aproximadamente al 1%, 5%, 10% y 20% del tamaño original:
par(mfrow=c(2,2))
plot(resize(grayscale(im),28,36),main="1%")
plot(resize(grayscale(im),143,177),main="5%")
plot(resize(grayscale(im),283,354),main="10%")
plot(resize(grayscale(im),567,709),main="20%")
Leer y transformar todas las imágenes
La librería imager contiene una serie de funciones para leer las imágenes de un directorio y convertirlas luego en un dataframe. Para ilustrar este proceso con un poco más de detalles se usará una función, leer_y_trs_img(), construida para mostrar cada paso de la transformación:
leer_y_trs_img<-function(img_name,path=NULL,x_l=10,y_l=10){
require(imager)
img_nombre<-paste0(path,img_name) # completa el nombre de la imagen con la ruta
imagen<-load.image(img_nombre) # carga la imagen
img_gris<-grayscale(imagen) # convierte la imagen a escala de grises
img_escalada<-resize(img_gris,x_l,y_l) # reescala la imagen
return(img_escalada)
}A continuación se muestra usa la función para cargar una sola imagen:
x<-leer_y_trs_img(filenames[1],path="./img/",x_l = 57,y_l = 71)
plot(x)
Ahora se procede a cargar todas las imágenes que se indican en filenames usando lapply():
lista_imagenes = lapply(filenames, leer_y_trs_img,path="./img/",x_l = 57,y_l = 71) El objeto lista_imagenes es una lista de imágenes. Cada componente es una imagen y se puede utilizar así:
par(mfrow=c(1,2))
plot(lista_imagenes[[2]])
plot(lista_imagenes[[8]])
Ahora se vectorizan las imágenes. Es decir, que vista la imagen como una matriz, sus columnas se ponen una debajo de la otra hasta obtener un vector. Estos vectores luego serán las filas de una matriz de datos donde cada pixel representa una variable y cada imagen representa una observación. El resultado de hacer esto para todas las imágenes es una matriz que tiene tantas filas como imágenes y tantas columnas como pixeles tengan las imágenes. La función as.numeric() aplicada sobre cada imagen devuelve un vector. Se aplica entonces la función as.numeric() sobre cada entrada del objeto lista_imagenes con la función lapply().
imagenes_vectorizadas<-lapply(lista_imagenes, as.numeric)
length(imagenes_vectorizadas[[1]]) # longitud del vector que representa la primera imagen## [1] 4047length(imagenes_vectorizadas) # cantidad de imágenes## [1] 72El resultado, imagenes_vectorizadas, es una lista de vectores que se concatenan llamando rbind() con do.call():
matriz_imagenes = do.call('rbind', imagenes_vectorizadas)
dim(matriz_imagenes) # dimensión de la matriz resultante## [1] 72 4047# as.data.frame.imlist es una alternativaEstadísticos descriptivos
Al tener una matriz de datos se pueden hacer todo lo que se hace con una matriz de datos. Sin embargo, el gran número de variables hace que utilizar la función summary() o la función pairs() sea poco práctico. Una aproximación puede ser generar estadísticas de resumen para cada pixel y luego ver estas estadísticas en forma de imágenes.
Veamos por ejemplo la imagen media:
imagen_media_vec<-apply(matriz_imagenes,2,mean) # calcula el valor promedio para cada pixel (columna)
length(imagen_media_vec) # longitud de la imagen promedio como vector## [1] 4047Ahora convirtamos la imagen promedio vectorizada en una imagen usando la función as.cimg():
imagen_media<-as.cimg(array(imagen_media_vec,dim=c(57,71)))
plot(imagen_media)
De igual manera se puede calcular la imagen que representa las desviaciones estándar para cada pixel:
imagen_sd_vec<-apply(matriz_imagenes,2,sd)
imagen_sd<-as.cimg(array(imagen_sd_vec,dim=c(57,71)))
plot(imagen_sd)
Descartando información poco relevante
La imagen de desviaciones estándares nos muestra varias zonas negras, donde el valor es cero. Estas zonas de baja variabilidad corresponden a pixeles poco informativos en esta muestra de imágenes. El fondo de la imagen, común para todos los individuos es, por ejemplo, una zona de baja variabilidad. Usemos el histograma de la imagen de desviaciones estándares para detectar estas zonas de baja variabilidad:
hist(imagen_sd)
Consideraremos que pixeles con una variabilidad inferior a 0.05 no serán relevantes en los siguientes análisis. La función threshold() asigna un 0 a cada pixel cuya intensidad está por debajo de cierto valor y 1 al resto de pixeles y nos devuelve una imagen. A continuación se usa esta función especificando un valor de umbral de 0.05:
imagen_sd_tr<-threshold(imagen_sd,thr = 0.05)
plot(imagen_sd_tr)
La imagen imagen_sd_tr es una máscara. Las zonas en blanco representan los pixeles que se incluirán en los análisis y las zonas en negro los pixeles que se descartarán por su baja variabilidad.
Veamos cuantos pixeles se descartan:
table(imagen_sd_tr)## imagen_sd_tr
## FALSE TRUE
## 739 3308Calculemos el porcentaje de información descartada:
table(imagen_sd_tr)[1]/sum(table(imagen_sd_tr))*100## FALSE
## 18.26044Esto quiere decir que aproximadamente el 18% de los datos están en zonas de baja variabilidad. Ahora identificamos los pixeles por fuera de las zonas de baja variabilidad usando la función which() y almacenamos la ubicación de estos pixeles en el objeto mascara:
mascara<-which(imagen_sd_tr==1)Ahora se extraen los datos correspondientes a los pixeles por fuera de la zona de baja variabilidad usando el objeto mascara para identificar las columnas relevantes en la matriz_imagenes. El resultado se almacena en el dataframe datos_con_mascara:
datos_con_mascara<-matriz_imagenes[,mascara]
dim(datos_con_mascara)## [1] 72 3308Análisis de componentes principales
A continuación se procederá a hacer una análisis de componentes principales sobre las imágenes. El primer paso será centrar y escalar la matriz de datos con la función scale(). Esto quiere decir que a cada columna se le resta su media y se divide por su desviación estándar. El resultado se almacena en el objeto datos_mascara_centrados:
datos_mascara_centrados<-scale(datos_con_mascara,center = TRUE,scale = TRUE)
dim(datos_mascara_centrados)## [1] 72 3308Notemos que nuestra matriz de datos tiene muchas más variables (3308) que observaciones (72). Esto hace que no se pueda usar la función princomp() para hacer el análisis de componentes principales. En su lugar se usa la función prcomp(), que utiliza la descomposición en valores singulares. A lo sumo se tendrán tantas componentes como observaciones (o filas de la matriz de datos):
modelo_pca<-prcomp(datos_mascara_centrados)
summary(modelo_pca)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 26.5193 25.1164 16.40629 15.33239 12.44363 9.93699
## Proportion of Variance 0.2126 0.1907 0.08137 0.07106 0.04681 0.02985
## Cumulative Proportion 0.2126 0.4033 0.48467 0.55573 0.60254 0.63239
## PC7 PC8 PC9 PC10 PC11 PC12 PC13
## Standard deviation 8.96227 8.51391 8.3537 7.88131 7.3197 6.72317 6.4567
## Proportion of Variance 0.02428 0.02191 0.0211 0.01878 0.0162 0.01366 0.0126
## Cumulative Proportion 0.65667 0.67858 0.6997 0.71846 0.7347 0.74832 0.7609
## PC14 PC15 PC16 PC17 PC18 PC19 PC20
## Standard deviation 6.11721 5.98341 5.6654 5.55175 5.29583 5.22120 4.97269
## Proportion of Variance 0.01131 0.01082 0.0097 0.00932 0.00848 0.00824 0.00748
## Cumulative Proportion 0.77223 0.78305 0.7928 0.80207 0.81055 0.81879 0.82627
## PC21 PC22 PC23 PC24 PC25 PC26 PC27
## Standard deviation 4.87272 4.76576 4.69384 4.58053 4.47333 4.32228 4.26341
## Proportion of Variance 0.00718 0.00687 0.00666 0.00634 0.00605 0.00565 0.00549
## Cumulative Proportion 0.83345 0.84031 0.84697 0.85331 0.85936 0.86501 0.87051
## PC28 PC29 PC30 PC31 PC32 PC33 PC34
## Standard deviation 4.20166 4.13948 4.09074 4.01823 3.93728 3.91091 3.74792
## Proportion of Variance 0.00534 0.00518 0.00506 0.00488 0.00469 0.00462 0.00425
## Cumulative Proportion 0.87584 0.88102 0.88608 0.89096 0.89565 0.90027 0.90452
## PC35 PC36 PC37 PC38 PC39 PC40 PC41
## Standard deviation 3.71081 3.67334 3.66160 3.55966 3.47681 3.46320 3.3525
## Proportion of Variance 0.00416 0.00408 0.00405 0.00383 0.00365 0.00363 0.0034
## Cumulative Proportion 0.90868 0.91276 0.91681 0.92064 0.92430 0.92792 0.9313
## PC42 PC43 PC44 PC45 PC46 PC47 PC48
## Standard deviation 3.32706 3.29960 3.18259 3.14377 3.10711 3.06952 3.06146
## Proportion of Variance 0.00335 0.00329 0.00306 0.00299 0.00292 0.00285 0.00283
## Cumulative Proportion 0.93467 0.93796 0.94102 0.94401 0.94693 0.94978 0.95261
## PC49 PC50 PC51 PC52 PC53 PC54 PC55
## Standard deviation 3.00845 2.95087 2.9350 2.89881 2.88403 2.84588 2.82912
## Proportion of Variance 0.00274 0.00263 0.0026 0.00254 0.00251 0.00245 0.00242
## Cumulative Proportion 0.95534 0.95798 0.9606 0.96312 0.96564 0.96808 0.97050
## PC56 PC57 PC58 PC59 PC60 PC61 PC62
## Standard deviation 2.76507 2.74264 2.69070 2.65895 2.60921 2.60380 2.52446
## Proportion of Variance 0.00231 0.00227 0.00219 0.00214 0.00206 0.00205 0.00193
## Cumulative Proportion 0.97281 0.97509 0.97728 0.97941 0.98147 0.98352 0.98545
## PC63 PC64 PC65 PC66 PC67 PC68 PC69
## Standard deviation 2.51638 2.48671 2.43214 2.38532 2.33871 2.28619 2.24074
## Proportion of Variance 0.00191 0.00187 0.00179 0.00172 0.00165 0.00158 0.00152
## Cumulative Proportion 0.98736 0.98923 0.99102 0.99274 0.99439 0.99597 0.99749
## PC70 PC71 PC72
## Standard deviation 2.14099 1.92727 1.178e-14
## Proportion of Variance 0.00139 0.00112 0.000e+00
## Cumulative Proportion 0.99888 1.00000 1.000e+00Las primeras 33 componentes principales explican aproximadamente el 90% de la variabilidad de los datos. Tomaremos los primeros 33 vectores propios (o en este caso vectores singulares) para formar las 33 primeras componentes principales:
eigen_modes_prin<-modelo_pca$rotation[,1:33]
dim(eigen_modes_prin)## [1] 3308 33Ahora se multiplicará cada imagen por los 33 vectores propios para representarla usando solo 33 variables:
datos_red<-datos_mascara_centrados%*%eigen_modes_prin
dim(datos_red)## [1] 72 33El resultado es una matriz de 72 filas y 33 columnas. Cada fila corresponde a un individuo y cada columna a una de las 33 primeras componentes principales.
Análisis de características en el espacio de componentes principales
Trataremos de ver si las en el espacio de las componentes principales es posible observar agrupaciones de imágenes por características. En la base de imágenes los individuos pueden ser jóvenes, adultos o adultos mayores, hombres o mujeres o tener una expresión facial específica: neutral (n), disgusto (d), miedo (f), felicidad (h), tristeza (s) o rabia (a).
Esta información esta codificada en el nombre de la imagen, veamos cómo:
head(filenames)## [1] "004_o_m_a_a.jpg" "004_o_m_a_b.jpg" "004_o_m_d_a.jpg" "004_o_m_d_b.jpg"
## [5] "004_o_m_f_a.jpg" "004_o_m_f_b.jpg"Así, las características pueden extraerse usando la función substr():
edad<-substr(filenames,start=5,stop=5)
sexo<-substr(filenames,start=7,stop=7)
expresion<-substr(filenames,start=9,stop=9)Veamos ahora cómo se agrupan las características en el espacio de las primeras cinco componentes principales usando la función pairs():
pairs(datos_red[,1:5],col=as.factor(edad), main="Edad")
pairs(datos_red[,1:5],col=as.factor(sexo), main="Sexo")
pairs(datos_red[,1:5],col=as.factor(expresion), main="Expresión")
¿Qué significa cada modo?
Veamos ahora la representación de cada vector propio. Para ello convertiremos cada vector propio en una imagen. Primero se crea una matriz que tendrá 33 filas, una para cada vector propio y \(57\times 71\) columnas. Usando la máscara definida anteriormente se llenan los pixeles (columnas) que sí tienen variabilidad:
modas<-matrix(0,ncol=57*71,nrow=33)
modas[,mascara]<-t(eigen_modes_prin)La función sweep() permite reescalar las modas multiplicándolas por las desviaciones estándar y luego sumándole la imagen media:
modas_escaladas<-sweep(modas,2,FUN='*',imagen_sd_vec)
modas_centradas<-sweep(modas_escaladas,2,FUN='+',imagen_media_vec)Ahora se convierte cada moda en una imagen:
modas_img_escalanat<-list()
modas_img<-list()
for (i in 1:33){
modas_img[[i]]<-as.cimg(modas[i,],x=57,y=71)
modas_img_escalanat[[i]]<-as.cimg(modas_centradas[i,],x=57,y=71)
}Finalmente visualizamos las primeras seis modas (sin reescalar):
par(mfrow=c(3,2))
lapply(X=1:6,FUN=function(x){plot(modas_img[[x]],main=paste0("Moda ",x))})
## [[1]]
## Image. Width: 57 pix Height: 71 pix Depth: 1 Colour channels: 1
##
## [[2]]
## Image. Width: 57 pix Height: 71 pix Depth: 1 Colour channels: 1
##
## [[3]]
## Image. Width: 57 pix Height: 71 pix Depth: 1 Colour channels: 1
##
## [[4]]
## Image. Width: 57 pix Height: 71 pix Depth: 1 Colour channels: 1
##
## [[5]]
## Image. Width: 57 pix Height: 71 pix Depth: 1 Colour channels: 1
##
## [[6]]
## Image. Width: 57 pix Height: 71 pix Depth: 1 Colour channels: 1Actividad
Construya un clasificador para las expresiones faciales usando regresión logística y análisis de componentes principales.